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Anharmonic phonons in few layer M0S2: Raman spectroscopy of ultra low energy 

compression and shear modes. 
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Molybdenum disulfide (M0S2) is a promising material for making two-dimensional crystals and 
flexible electronic and optoelectronic devices at the nanoscale^^. M0S2 flakes can show high mobil- 
ities and have even been integrated in nanocircuits^^. A fundamental requirement for such use is 
efficient thermal transport. Electronic transport generates heat which needs to be evacuated, more 
crucially so in nanostructures. Anharmonic phonon-phonon scattering is the dominant intrinsic 
limitation to thermal transport in insulators. Here, using appropriate samples, ultra-low energy 
Raman spectroscopy and first principles calculations, we provide a full experimental and theoretical 
description of compression and shear modes of few- layer (FL) M0S2. We demonstrate that the com- 
pression modes are strongly anharmonic with a marked enhancement of phonon-phonon scattering 
as the number of layers is reduced, most likely a general feature of nanolayered materials with weak 
inter layer coupling. 
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I. INTRODUCTION 

Bulk M0S2 is made of vertically stacked layers (single 
formula unit consisting of a Mo sheet sandwiched be- 
tween two S sheets) weakly held together by Van der 
Waals forces, with 2 layers per unit cell. While M0S2 
is an indirect gap semiconductor, it displays a crossover 
to a direct gap semiconductor with a resulting marked 
increase of photoluminescence upon reduction of crystal 
thickness down to one layei^. 

In the case of a layered hexagonal system, the shearing 
modes are twofold degenerate as there are two equiva- 
lent in-plane shear-directions. Generally an N-layer flake 
has 2 * {N — 1) shearing modes and N — 1 compression- 
extension (noted compression henceforth) modes. Thus, 
while single layer M0S2 should have neither one nor the 
other mode, in FL M0S2 additional modes should ap- 
pear with respect to bulk or the bilayer which both have 
2 degenerate shear modes and one compression mode. 

Both modes, importantly, are expected at very low en- 
ergies since the interlayer interaction in M0S2 is weak 
{uUiy < 55 cm~^ or 74K). Therein lies their relevance to 
transport, both thermal and electronic. At room temper- 
ature these phonon modes are all thermally populated 
and influence thermal transport via phonon scattering 
with defects and impurities or via phonon-phonon scat- 
tering. As anharmonicity is an intrinsic mechanism, in 
clean samples it is the dominant limitation to the thermal 
conductivity. Determining the behaviour of these modes 
as a function of the flake thickness is of the greatest im- 
portance for nanoelectronic devices based on M0S2. 

Recently low energy Raman modes in Graphene^ and 
MoSj^H^were measured for FL flakes on p-doped Si sub- 
strates. In Ref. '9' only a single shear mode was detected 
and no compression modes were seen. Low energy modes 
are weak in intensity and, to eliminate a broad low-energy 
background due to inelastic scattering of free carriers in 




FIG. 1: The few-layers M0S2 flake made by anodic bonding 
and used for the experiment. The number of layers is indi- 
cated. 



the p-doped substrate^, these spectra were recorded in 
crossed geometry (i.e. the polarization of the outcom- 
ing light is perpendicular to that of the incident beam) 
where all compression modes are forbidden. In Ref. [TTJ 
beside the shearing mode detected in Ref. 9, a second 
feature was detected. No polarization analysis of the Ra- 
man spectra was performed and one of the features was 
attributed to a compression mode from a flt to a 1/N 
behaviour {N = number of layers). In Ref. [10] a detailed 
study of shearing and compression modes up to 19 layers 
has been carried out. The authors were able to classify 
the low energy Raman peaks in two groups, namely those 
that stiffen with increasing N and those that soften with 
increasing N. Then using a chain model they were able to 
build fan diagrams and obtain M0S2 shearing and com- 
pression strength. No flrst principles calculations of the 
Raman spectra were carried out. 
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FIG. 2: Energy difference of the E2g and the Ai^ phonon modes at high energy versus the number of layers (left) and 
experimental high-energy Raman spectra as a function of the layer number (right). 



In this work we measure low-energy Raman spectra as 
a function of the number of layers in multilayer M0S2 • We 
overcome the difficulty related to the presence of a broad 
low-energy background due to inelastic scattering of free 
carriers in the p-doped substrate^ by using FL M0S2 on 
borosilicate glass substrates. We measure shear and com- 
pression modes from 1 to 5 layers. By performing first 
principles calculations of the position^ and intensity^ 
of Raman peaks we obtain a complete understanding of 
shear and compression modes in FL M0S2 • We also an- 
alyze theoretically the dependence of the main shearing 
mode as a function of applied pressure and show that it 
behaves linearly at low pressure (below IGPa). 



II. EXPERIMENTAL 

The samples are made using the anodic bonding 
method which bonds a bulk flake to a borosilicate glass 
substratel^nS ^t a temeprtaure between 130C and 200C 
and a high voltage which may range from 200V to 1500V. 
The flake is then mechanically cleaved leaving large few 
layer samples on the glass (see Fig. IT]). The sample used 
in this work is exceptional in that it provides all different 
thicknesses used on the same flake, making comparison 
easy. The sample thickness was identified first by optical 
contrast and then confirmed by atomic force microscope 
and Raman spectroscop}'^^, as shown in Fig. [2] 

MicroRaman spectra of the FL M0S2 are measured 
(532 nm radiation, ^ 1 mW laser power) in backscat- 
tering configuration with parallel and crossed polariza- 
tion geometry. To enable measurements down to ~ 10 
cm~^ on a single-grating spectrometer (Lab RAM HR 
from HORIBA Jobin Yvon), an ultra- low wavenumber 
filtering (ULF-^^) accessory for 532 nm wavelength was 
used. These volume Bragg gratings can be fabricated 



with diffraction efficiencies as high as 99.99% and the 
linewidth narrower than 1 cm~^ at FWHM that corre- 
sponds to 3-4 cm~^ cut-off frequency at -60 dB from max- 
imum. They also make a unique notch filter for Rayleigh 
light rejection by sequential cascading of several Bragg 
notch filters, enabling ultra- low frequenc y Ram an mea- 
surements with single stage spectrometers^^Ell Xo avoid 
laser heating, a laser power of 1 mW was focused through 
a 50X or lOOX microscope objectives. 

In order to evaluate the actual temperature on the 
sample and to exclude laser-heating, we measure for all 
thicknesses both Stokes and antiStokes features at high 
and low energy. We then determine the temperature as 
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where u is the Raman shift and ul the pulsation of the 
laser light. For a more precise determination of the tem- 
perature we use the intensity of the high energy Ai^ and 
E2^ modes. We obtain T ^ 360K. Thus the laser heating 
is ^ 60 K, i. e. negligible. 



III. THEORY 

Calculations were performed by using density func- 
tional theory in the local density approximation^^. The 
QUANTUM-ESPRESSO^^ package was used with norm- 
conserving pseudopotentials and a plane- wave cutoff en- 
ergy of 90 Ry. Semicore states were included in Mo pseu- 
dopotential. All calculations were performed at and 
6 Kbar uniaxial pressure, corresponding to the uniaxial 
pressure imposed in the anodic bonding procedure. The 
crystal structure at a given pressure is obtained by keep- 
ing the same in-plane experimental lattice parameter as 
in bulk samples. The interlayer distance is obtained by 
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FIG. 3: Experimental (blue) and theoretical (red-dashed) Raman spectra in parallel (left) and crossed (right) configuration. 
The inset shows a blow up of the low energy region. Also shown are the schematized M0S2 layer movements for the compression 
(left) and the shear (right) modes. All arrows may be simultaneously reversed and the size corresponds to the amplitude 



imposing a 6 Kbar pressure with respect to the bulk ex- 
perimental structure at Kbar. The distance between 
the plane is then kept constant for all N-layers flakes, but 
Sulfur height, the only free parameter, is optimized fol- 
lowing internal forces. Phonon frequencies, born-effective 
charges and Raman tensor were calculated using a 8 x 8 x 1 
k-point grid for the monolayer and a 8 x 8 x 8 for the 
bulk. Using this grid, phonon frequencies of shearing 
and compressing modes are converged with an accuracy 
of 1 cm~^. Raman intensities were calculated with the 
method of Ref. 21, in the Placzek approximation. 

The intensity of a mode u is written as 

r oc I^ (n, + l)/cj, 

where oOi, and riu are the phonon frequency and the occu- 
pation of the phonon mode u. Moreover, /q = |e^ • Ae^p 
where A is the Raman tensor while e^ and e^ are the the 
polarization of the incident and scattered radiation re- 



The experimental spectrum is then obtained as 



spectively. In table III we give the calculated value of Iq 
as a function of the layer number and the experimental 
geometries for the low energy modes. 
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In order to compare with experiments, the Dirac S func- 
tions are convoluted with the experimental linewidths. 

At ambient pressure, the calculated frequencies are in 
excellent agreement with previous calculations!^, how- 
ever they disagree with the calculations of Ref. Hi] In 
particular, for a M0S2 bilayer at zero uniaxial pressure 
and using the experimental in-plane lattice parameter we 
find 20 cm~^ for the main shearing mode. Performing 
structural optimization of both in-plane and interlayer 
distance (keeping always the same empty region between 
periodic images) we find 24.49 cm~^ for the shear mode. 
Thus the shear phonon frequency weakly depends on the 
choice of the experimental or theoretical in-plane lattice 
parameters in the calculation. This has to be compared 
with 22 cm~^ in our experimental Raman data, with ^ 23 
cm~^ in experimental Raman data of Ref. [TTIand with 
^ 19.5 cm-^ in Raman data of Ref. H In Ref. ^ 
the shearing mode was calculated at 35.3 cm~^ using the 
theoretical lattice structure and the LDA approximation 
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TABLE I: Raman intensity Iq for different modes ly and as a 
function of the number of layers N for backscattering geome- 
try (labeled || geom.) and cross cackscattering geometry. The 
intensity is normalized to the most intense low energy mode 
(below 100 cm~^ ). Note that in the two and three layer case, 
the compression mode has stronger intensity Iq then the main 
shear mode, however the larger linewidth suppresses its /^. 



innermost layers, these being more tightly bound by the 
bilateral interaction with the other layers. The hardest 
(also named principal or primary) shear mode is labeled 
Sn, where n is the number of layers. The shear mode Sn 
is Raman active and is also the most intense of all the 
shear modes. As the number of layer increases, there are 
more possibilities of rigidly shifting layers. For example 
in a three layer flake, it is possible to shift only the top 
(or bottom) layer keeping the other two fixed. This mode 
is however softer then the main shear mode, as the outer 
layers have only one nearest neighbour layer with weaker 
binding. In our case and at 6 Kbar uniaxial pressure, the 
secondary shear mode is calculated to be at 16.3 cm~l, 
but with essentially zero Raman intensity (although this 
mode is not forbidden by symmetry). In the general case 
of an N-layer flake, there are {N — 1) independent ways of 
shifting a subset of layers with respect to the others and 
their energies lie between that of the (softest) secondary 
shear mode related to the shift of extremal layer and the 
(hardest) main shear mode. 

Compression (or extension) modes are rigid vibrations 
of the layers in the direction perpendicular to the layers. 
For an N layer flake (with N< 5) the rigid layer displace- 
ment pattern are schematically illustrated in Fig. [3] (left). 
As in the case of shear modes it is possible to identify a 
primary or main compression mode. The main compres- 
sion mode is the one corresponding to (i) the lower half 
of the layers shifting in the same direction and (ii) the 
higher half of the layers in the opposite direction. The 
main compression mode (labeled Cn) is the softest com- 
pression mode. Secondary compression modes arise when 
the top half (respectively bottom half) layers are not all 
displaced in the same direction. Secondary compression 
modes are higher in energy then the main compression 
mode (see Fig. |4J the rigid layer model of ^23^ and the 
discussion below^ 
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which would correspond to a very large applied pressure, 
not relevant to the experiments in consideration. 



IV. GENERAL DISCUSSION ON SECONDARY 
SHEAR AND COMPRESSION MODES 

As already stressed in the introduction, a M0S2 N-layer 
fiake has 2 * {N — 1) twofold degenerate shearing modes 
and N — 1 compression modes. Shear modes correspond 
to rigid layer displacements perpendicular to the c-axis. 
The twofold degeneracy depends on the crystal symmetry 
of the lattice. In M0S2 for example it is equivalent to 
rigidly shift a subset of layers with respect to one in- 
plane crystalline axis or the other. For an N layer flake 
(with N< 5) the rigid layer displacement patterns are 
schematically illustrated in Fig. [s] (right). The hardest 
shear mode always corresponds to the rigid shift of the 



In Fig. [3] we show measured and calculated spectra in 
parallel and crossed geometries. The peaks are normal- 
ized to the main shear mode intensity and the theoretical 
spectra are convoluted by the experimental linewidth. 

By comparing the main shear- mode (labeled Sn) to ex- 
isting measurements^, we notice that in our samples this 
mode is systematically harder by ^ 3 cm~^. We attribute 
this to the anodic bonding method which binds flakes 
to a glass substrate electrostatically due to the creation 
of a space charge in the substrate. This also generates 
an electrostatic pressure on the flake bound to the sub- 
strate. Knowing the depth of this space charge layer (1-2 
microns) this uniaxial pressure along the c-axis^^ can be 
estimated to be in the range of 3-6 Kbar. 

We thus performed first principles calculations as a 
function of pressure and obtain an essentially linear be- 
haviour of the main shearing mode for a M0S2 bilayer as a 
function of uniaxial pressure, with a linear coefficient that 
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FIG. 4: Fan diagram of calculated shear (red) and compres- 
sion (black) mode frequencies as a function of the number of 
layers. Full symbols represent modes visibles in our Raman 
experiments in parallel and crossed configuration. The labels 
of the different modes are the same as in Fig. [3] 



is 1.17 cm~^/Kbar (the second order coefficient in the fit 
is 0.04 cm~^/Kbar^). Comparing the measured spectra 
with the calculated ones for 6 Kbar uniaxial pressure, we 
find remarkable agreement. 

Besides the main shear mode, already detected in Refs. 
19)11! we measure secondary shear (Sn and S^) and com- 
pression (labeled C^ and C^) modes (see Figs. [4]and|3]). 
In Ref. [11] only the Cn compression mode was detected. 

The compression mode is clearly visible in 2, 4 and 5 
layer samples and less so in the 3 layer sample because 
it coincides in energy with the shear mode. In the 4 and 
5 layer samples, theory also accounts for the additional 
shear (S^) and compression (C^) modes detected in ex- 
periments. 

It is worthwhile to recall that the energy of shear and 
compression modes can either increase or decrease as a 
function of the number of layers A/", as shown in the fan 
diagrams in Ref. 23 and in Fig. gj In FL MoS2, the 
energy of Raman visible E^^ shear modes increases with 
the number of layer, while that of Raman visible A\g 
compression modes decreases. 

The qualitative behaviour of the position of the main 
shear and compression modes as a function of layer num- 
ber can be easily understood in a simple rigid-layer (or 
chain) modeP^. We write for the Sn shear- mode fre- 
quency 
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V 27rc V M 
where \i = 30.75 Kg/m^ is the rigid layer mass per unit 
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of compression (top) and shear (bottom) modes as a function 
of the inverse layer number. 



and OLs is the shearing strength. The Cn compression 
mode behaves as 
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and OLc is the compression strength. From the experi- 
mental curves we ob tain a 5 = 27.44 x lO^^N/m^, twice 
the value in graphen^^l^, and ac = 420.44 x lO^^N/m^. 
Both these effects are in part explained by the smaller 
interlayer distance in M0S2 (the Sulfur-Sulfur distance 
along c is 3.03A) with respect to graphite (3. 35 A). The 
agreement between Eqs. [3J |4] and experimental data is 
shown in Fig. |5J validating the chain model and the ex- 
tracted values of as and ac- 

In Fig. |5] we also plot the variation of the linewidth 
7 for the Sn shear and the Cn compression modes as a 
function of the number of layers n. Shear mode linewidth 
is resolution limited while all the compression modes are 
very broad (roughly 7 times broader), the linewidth of 
the C2 mode being the largest. As the broadening is 
inversely proportional to the phonon-phonon scattering 
time, our result indicates that the phonon scattering time 
of compression modes is approximately 7 times smaller 
than that of shear modes. Thus the contribution of opti- 
cal modes to the intrinsic thermal conductivity of M0S2 
flakes is dominated by scattering to compression modes. 



In a M0S2 bilayer, the scattering time of compression 
modes is nearly 9 times smaller then that of shear modes. 
The compression modes linewidths are generally larger 
then those of shearing modes as the potential is more an- 
harmonic for a displacement perpendicular to the M0S2 
layers then for a shearing displacement. In the case of C2, 
the linewidth is enhanced with respect to C^ with n > 2 
due to the fact that more channels for anharmonic decay 
are available. Indeed for n > 2, the compression mode is 
lower in energy (or at the same energy for n = 3) then the 
shearing mode. As such it can only decay into two acous- 
tic modes of opposite momentum. In the case n = 2, the 
compression mode is at roughly twice the energy of the 
shearing mode 82- Thus the compression mode can decay 
into (i) two acoustic modes of opposite momentum, (ii) 
an acoustic and a shearing mode of opposite momentum, 
or (iii) two shearing modes of opposite momentum. 



strongly anharmonic, with phonon-phonon scattering in- 
creasing upon reducing the number of layers. Thus com- 
pression modes represent the overriding optical phonon 
contribution to the intrinsic thermal conductivity of 
M0S2 flakes, a crucial aspect of any use of these in future 
nano or microelectronic devices. The relevance of our 
work is far reaching as compression modes are most likely 
strongly anharmonic in all flakes obtained from weakly- 
interacting layered- materials such as few layer graphene, 
transition metal dichalcogenides and topological insula- 
tors. In all these systems a crucial limit to thermal trans- 
port could be the anharmonicity of compression phonon 
modes. 
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VI. CONCLUSION 

In conclusion we have measured primary and sec- 
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in concmsion we nave measurea primai 
ondary shearing and compression modes in M0S2 from 
1 to 5 layers. The compression modes are found to be 
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